Poloidal-toroidal decomposition in a finite 
cylinder. I. Influence matrices for the 
magnetohydrodynamic equations 



Piotr Boronski and Laurette S. Tuckerman 

LIMSI-CNRS, BP 133, 91403 Orsay, France 



Abstract 

The Navier-Stokes equations and magnetohydrodynamics equations are written in 
terms of poloidal and toroidal potentials in a finite cylinder. This formulation in- 
sures that the velocity and magnetic fields are divergence-free by construction, but 
leads to systems of partial differential equations of higher order, whose bound- 
ary conditions are coupled. The influence matrix technique is used to transform 
these systems into decoupled parabolic and elliptic problems. The magnetic field 
in the induction equation is matched to that in an exterior vacuum by means of the 
Dirichlet-to-Neumann mapping, thus eliminating the need to discretize the exte- 
rior. The influence matrix is scaled in order to attain an acceptable condition num- 
ber. 



1 Motivation and Governing Equations 

The requirement that velocity and magnetic fields be solenoidal, i.e. divergence-free, 
represents one of the most challenging difficulties in hydrodynamics and in magneto- 
hydrodynamics [1, 2, 3, 4, 5, 6, 7]. For the velocity field, this condition is the fundamen- 
tal approximation used in incompressible fluid dynamics. For the magnetic field, this 
condition is the statement of the non-existence of magnetic monopoles. 

Two main approaches exist for imposing this requirement. The first is to use three field 
components and to project three-dimensional fields onto a divergence-free field. In an 
incompressible fluid, the pressure serves to counterbalance the nonlinear term which 
is the source of the divergence in the Navier-Stokes equations; the pressure also plays 
this role numerically. The divergence of the Navier-Stokes equations is taken, lead- 
ing to a Poisson problem for the pressure. However, the boundary conditions on the 
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equations for (u,p) involve only the velocity, leading to coupling between the equa- 
tions to be solved for u and p [3, 4]. The coupled equations can be solved in several 
stages by a Green's function or influence matrix method [3]. In projection-diffusion 
schemes, approximate boundary conditions are imposed for the pressure [5]. For mag- 
netic fields, however, the exact evolution of the equations conserves divergence and 
there exists no analogue to the pressure. Thus if the numerical algorithm creates di- 
vergence, there is no mechanism for eliminating it and it may accumulate [6]. For this 
reason, magnetohydrodynamic codes sometimes include a fictitious term analogous to 
the hydrodynamic pressure, which must be treated numerically [7]. 

The second approach, which is the focus of this paper, is to express fields in such a way 
that they are divergence-free by construction. It can be proved that a field F which is 
solenoidal (divergence-free) in a simply connected domain can be written as: 

F = V x (ipe) + V x V x (cpe) (1.1) 

where e denotes a unit vector. In addition to being divergence-free, F has the advan- 
tage of involving only two scalar fields. This makes more economical use of computer 
memory and allows all calculations to be implemented using only scalar fields. 

Equations governing the evolution of the two potentials are derived by taking the curl 
and double curl of the original equations, increasing the order of the differential equa- 
tions. In addition, boundary conditions, some also of high order, couple the two poten- 
tials. In certain geometries with two periodic directions, these are only minor obstacles 
[1]. In spectral treatments of such geometries, the basis functions insure periodicity, 
which is preserved under differentiation and addition. At most, special consideration 
must be given to constant modes. The standard examples are a spherical geometry 
[8, 9, 10, 11, 12, 13] or a three-dimensional Cartesian geometry with one bounded di- 
rection and two perpendicular periodic directions, such as channel flow [14, 15]. Other 
applications are in a cylindrical geometry with periodic z and 8 directions [1, 16, 17]. 

In geometries with more than one nonperiodic direction, far more care is required. 
Marques [1] gave a detailed analysis of the poloidal-toroidal decomposition for the 
Navier-Stokes equations and its formulation and validity for general topologies. This 
analysis was then put into practice in a linear stability analysis of Rayleigh-Benard 
convection in a finite cylindrical geometry [2]. However, the governing equations de- 
rived in [2] contain large linear systems that couple the potentials and their laplacians 
and bilaplacians, but whose solution would be required in implicit time integration. 
Analogous problems arise in the other formulations of incompressible fluid dynam- 
ics. In the 2D streamfunction-vorticity formulation, the equations for the vorticity and 
the streamfunction are coupled by the fact that boundary conditions exist only for the 
streamfunction and none on the vorticity. In the (u,p) primitive variable formulation, 
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the pressure is the solution to a Poisson problem for which the appropriate boundary 
condition is that the velocity be divergence-free [3, 4]. 



Our primary goal in this paper is to demonstrate that the high-order equations can be 
separated via the influence matrix technique into a sequence of problems of lower or- 
der, each with its own boundary conditions, as was done for the primitive variable for- 
mulation in [3]. This makes implicit time integration feasible for the poloidal-toroidal 
decomposition in geometries with two non-periodic directions. A secondary goal is 
to carry out the same analysis for a magnetic field which is governed by the induc- 
tion equation and which generalizes the Navier-Stokes equation by the inclusion of 
the Lorentz force. 

The equations we will consider are the magnetohydrodynamic equations: 

g2 

9fii + (u • V)u = (B • V)B + Re _1 Au - V(p + y ) (1.2a) 

V • u = (1.2b) 

d t B = V x (u x B) + Rra" 1 AB (1.3a) 

V • B = (1.3b) 

where Re is the usual hydrodynamic Reynolds number and Rm the magnetic Reynolds 
number. Equations (1.2) and (1.3) are of different types: for a divergence-free magnetic 
field B, all the terms of (1.3) have zero divergence as well, but this is not the case for 
(1.2). 

The velocity and magnetic fields are to be calculated in a finite cylinder. We consider 
specifically the case in which the flow is driven by rotating upper and lower disks, 
although our method does not depend on this. For disks rotating in opposite directions 
this configuration is called the von Karman flow [18, 19, 20]. The magnetic field inside 
the cylinder is required to match the field outside, which goes to zero at infinity. These 
boundary conditions are expressed as: 

u = at r = 1, (1.4a) 

u = rco±eg at z = ±-, (l-4b) 

B int - B ext = on 90, (1.5a) 

B = at infinity. (1.5b) 

where O denotes the interior domain (the cylinder) and dCl is its boundary. 

The poloidal and toroidal components for this configuration in the axisymmetric case 
with e = e 2 are illustrated in figure 1. The toroidal flow corresponds to motion with 
only azimuthal velocity. The poloidal flow forms recirculation rolls in the (r,z) plane. 
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Fig. 1. Axisymmetric flow between counter-rotating disks. Poloidal component: solid curves. 
Toroidal component: dashed curves. 

For a non-axisymmetric flow, there is no clear correspondence between each potential 
and a simple topological structure. 

In section 2, we give a general description of the poloidal-toroidal decomposition. In 
section 3, we then specialize to the Navier-Stokes equations in a finite cylinder, formu- 
lating the boundary conditions for this case. In section 4 we show how to decouple the 
equations and boundary conditions via the influence matrix technique. Finally, in sec- 
tion 5, we present the equations and boundary conditions for the induction equation 
which governs the magnetic field, and the corresponding influence matrix. 



2 Poloidal-toroidal decomposition 



2.1 Governing equations 



The poloidal-toroidal decomposition generalizes to three dimensions the two-dimen- 
sional streamfunction-vorticity formulation. We follow the analysis and notation of 
[1], but specializing to the case of a domain which is contractible to a point (i.e. has no 
holes). Then: 

V • F = F = VxA (2.1) 

A distinguished direction and associated unit vector e is selected and A can be decom- 
posed such that: 

F = Vxfe + VxVxfe (2.2) 

The direction e is called vertical and those perpendicular to e are called horizontal; 
see figure 2. A number of possibilities exist for e. Among these, the choices e = e z (in 
Cartesian or cylindrical coordinates) or e = (the spherical radius) decouple xp and <p 
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in the diffusive operators since: 



e F = -A h (j), e • AF = -AA h (p, (2.3a) 

e • V x F = -A h xp, e • V x AF = -AA h xp, (2.3b) 

e • V x V x F = AA h <p, e • V x V x AF = AAA h <p. (2.3c) 

where A^ is the two-dimensional Laplacian acting in the horizontal directions, i.e., 
those perpendicular to e. (The decoupling (2.3) does not hold [1] when the cylindri- 
cal radius e r is chosen as the distinguished direction e). 

The equations for the velocity potentials are derived by taking the e component of the 
single and double curl of (1.2); those for the magnetic potentials are derived by tak- 
ing the e component itself and the single curl of (1.3). The difference arises from the 
fact that all the terms of (1.3) are divergence-free and there is no pressure to eliminate. 
Combining (1.2)-(1.3) and (2.3) leads to the evolution equations for the scalar poten- 
tials: 

(d t - Re- 1 A)A h ip u = e • V x S M (2.4a) 
(d t - Re~ 1 A)AA h (p u = -e • V x V x S M (2.4b) 

(d t - Rm- 1 A)A h <p B = e • S B (2.5a) 

(d t - Rm- 1 A)A h tp B = e • V x S B (2.5b) 

where: 

S M = (u • V)u - (B • V)B (2.6a) 

S B = - V x (u x B) (2.6b) 

Equations (2.4)-(2.5) are not all of the same order in the vertical and horizontal di- 
rections. For example, for the velocity, (2.4a) is 2 nd order in the vertical direction and 
4 th order in the horizontal directions, while (2.4b) is 4 th order in the vertical direction 
and 6 th order in the horizontal directions. A corresponding number of boundary con- 
ditions are required for the velocity potentials, a total of (2+4) /2=3 conditions at each 
vertical boundary and (4+6) /2=5 at each horizontal boundary for u. The conditions at 
the vertical boundaries are those corresponding to the physical problem. At the hor- 
izontal boundaries, the physical conditions must be supplemented by two additional 
conditions whose derivation is the subject of the remainder of this section. 

2.2 Gauge freedom 

The poloidal-toroidal formulation (2.2) contains a gauge freedom for the choice of ip 
and <p, which is identified by finding the class of potentials satisfying the homogeneous 
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Fig. 2. Geometry for potential variable formulation. Q is a cylindrical domain. The vector e 
points in the distinguished vertical direction, here e z . Q/, are slices of Q perpendicular to e, 
here disks. The boundary of Q/, is 30/,, here a circle. The vector n is normal to both e and to 
f)/,; here n = e r . 

problem F = 0. For e = e z (Cartesian or cylindrical coordinate) or e = e^ (spherical 
radius), this leads to: 



jhom 



= V x (V om e) +V x V x (V om e) 
= ex V h ip hom +V h d e (p hom - (kh<p hom ) 



e • Y hom = 



A h( p hom = 

horn 



e x F" om = => ex V h ip nom = -V h d e (p' 



Jiom 



(2.7a) 
(2.7b) 



where e denotes the coordinate corresponding to e. The existence of tp hom satisfying 
(2.7b) for all cp hom satisfying (2.7a) is demonstrated as follows. Condition (2.7a) implies: 



V h -(V h d e <p hom )=0 



(2.8) 



Applying (2.1) to the simply-connected two-dimensional domain slices perpendicular 
to e, (2.8) implies that there exists a ip hom satisfying: 



V h d e (p' 



horn 



V h x (-ip hom e) = e x V h ip 



horn 



(2.9) 



Thus, the poloidal potential (p is determined up to a harmonic function on each domain 
slice perpendicular to e: 



cp^(p + ( p h ° m ■ A h( p hom = 



(2.10a) 
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while xp is determined up to an arbitrary function of the coordinate e: 

ip~ip + h(e) (2.10b) 
The choice of gauge constitutes one of the two additional conditions required. 



2.3 Compatibility condition 



We have not yet demonstrated the equivalence between the potential and primitive 
variable formulations. Since the curl of equations (1.2) and (1.3) were taken, they gained 
an additional degree of freedom which must be fixed in such a way that these equa- 
tions in potential form (2.4)-(2.5) define the same velocity u and magnetic field B as the 
original MHD equations (1.2)-(1.3). We will require the fact that on a simply-connected 
domain, a field is a gradient if and only if it is curl-free: 

£ = Vp ^ Vxf = (2.11) 

which is a consequence of Stokes' theorem. We will first write (1.2)-(1.3) in a compact 
form, which will let us use a common form for (2.4) and (2.5): 

f u = (d t - Re- 1 A) u + S H = -V (p + B 2 /2) (2.12a) 
g M = Vxf„ =0 (2.12b) 

g B = (d t - Rm- l ^j B + S B = (2.12c) 

where (2.12a) and (2.12b) are equivalent, by (2.11). Then we can write the primitive 
variable formulation (1.2)-(1.3) and potential formulation (2.4)-(2.5) using for either 
g = g M org = g B : 

primitive variables potential formulation 

e • g = (2.13) 



g = 



e- V x g =0 



Marques [1] proves that in a simply connected domain O, the potential and primitive 
variable formulations are equivalent if additional conditions are satisfied: 



g = & 



e • g = in O (2.14a) 

e • V x g = in O (2.14b) 

V-g = inO (2.14c) 

n • g = on dC} h (2.14d) 
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In (2.14d), n is the vector normal to the boundary of slices perpendicular to e. We 
recall that in our case, e = e Z/ the slices Ci^ are disks, their boundaries 90^ are circles, 
and n = e r is the radial unit vector. We illustrate this geometry in figure 2. 

The rightwards implication of (2.14) is obvious. The leftwards implication of (2.14) is 
proved as follows. We first use (2.14a) and the two-dimensional Stokes' Theorem (2.11) 
to introduce a scalar function k 



g = Vfc/c (2.15) 



recalling that the subscript h restricts differential operators to the directions perpendic- 
ular to e; in our case the (r,9) directions. We then use the additional divergence-free 
condition (2.14c) to show that K is harmonic: 




A h K = (2.16) 



v-g = o 

The additional condition (2.14d) then provides a Neumann boundary condition on k: 



K = K (e) (2.17) 



A h K = 
n • Vk = on dCl^ 

(see figure 2, where e = z). Finally 

g = V h K(e)=0 (2.18) 

since the V/, measures variation in the horizontal directions, which are perpendicular 
to the coordinate e. The divergence-free condition (2.14c) is satisfied for u since (2.12b) 
defines g M as a curl. It is satisfied for B because (2.6) is divergence-free if B is, e.g. if B 
is expanded as (2.2). 

Condition (2.14d), which is called the compatibility condition and which ensures the 
equivalence of both formulations, is the projection of the original equations normal to 
the boundary. Its interpretation is quite intuitive: the compatibility condition preserves 
information about the original equations which has been lost by taking the curl. This 
procedure is familiar from simpler contexts: when an equation is differentiated, it must 
be supplemented by a constant of integration, which is the evaluation of the original 
equation at a point. Condition (2.14d) is sufficient but not unique - other boundary 
conditions ensuring (2.14) exist. 

We can extend the equivalence (2.14) proved in [1] to justify the transformed boundary 
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conditions often used in practice in the toroidal-poloidal formulation. To impose the 
boundary condition u — u bc = on a simply connected boundary dCl with normal e, 
we substitute for the three vector components the conditions: 

e • (u - u bc ) =0 on 90 (2.19a) 

e • V x (u - u bc ) =0 on 90 (2.19b) 

Vfe • (u - u bc ) =0 on an (2.19c) 

n (u - u bc ) =0 on 3(dfi) (2.19d) 

where 3 (an) is the one-dimensional boundary of an and n is perpendicular both to 
this boundary and to e. Equation (2.19c) can be replaced by 



V • [(e • u)e] = V • u - V fe • u bc = on an (2.19c') 

where the second equality is valid when u is divergence-free and the boundary condi- 
tions are homogeneous. 

The transformed boundary conditions (2.19) are familiar in the context of a spherical or 
infinite planar surface, where the additional condition (2.19d) is not needed since these 
surfaces have no boundaries. For example, in the case of flow between two stationary 
infinite planes at z = ±1, boundary conditions (2.19) take the form: 



w = atz = ±l (2.20a) 
tj = atz = ±l (2.20b) 
d z w = atz = ±l (2.20c) 

where w and r\ are the vertical velocity and vorticity. 

The derivation of both the gauge and the compatibility conditions depend on prop- 
erties of the horizontal Laplace equation; see (2.10a) and (2.17). If the only harmonic 
function is a constant, the Neumann condition in (2.17) is superfluous. This is the case, 
for example, on the surface of a sphere. All horizontal directions are periodic, so that 
functions which grow monotonically in these directions are excluded; the compatibil- 
ity condition (2.14d) can simply be dropped. However, in domains with more compli- 
cated topologies, such as those bounded by two infinite planes or cylinders considered 
to be doubly periodic, additional conditions are necessary in order for (2.11) to hold. A 
derivation of these conditions for a general domain can be found in [1]. 
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3 Conditions on the velocity field 



We now turn to the conditions to be imposed on the velocity field in the finite-cylindrical 
geometry for which e = e z ; see [2]. Because the next two sections will refer exclusively 
to the velocity we drop the subscript u. For reference, we write for F defined in (2.2) 
the identities: 

F = — e z x V h ip + V h d z (p - e z A h (p (3.1a) 
V x F = -e z x V h A(p + V h d z xp - e z A h xp (3.1b) 
AF = -e z x V h Atp + V h d z A(p - e z A h A(p (3.1c) 
V x AF = e 2 x V^AA^ + V h d z A\p - e z A h Axp (3. Id) 

which will facilitate calculations of vector quantities. 

3.1 Gauge and boundary conditions 

The governing equations are: 

(d t - Re~ 1 A)A h tp = e 2 • V x S (3.2a) 
(d t - Re~ 1 A)AA h (p = -e z • V x V x S (3.2b) 

The system (3.2) contains five Laplacians acting in the horizontal directions and three 
acting in the vertical directions. Three conditions in each direction are derived from 
the velocity boundary conditions. The two remaining conditions in the horizontal di- 
rection are the gauge and compatibility conditions. 

The simplest choice of gauge is: 

= at r = l (3.3a) 

along with 

^ = at r = (3.3b) 

On the cylinder, boundary conditions are imposed on u r , uq, u z . Referring to (3.1a), we 
have: 

1 

U r = -dglp + d rz (p = 
1 

Uq = -d r lp + -d 9z (p = 



r 

u z = -A h cp = 



(3.4a) 

at r = 1 



(3.4b) 
(3.4c) 



10 



The gauge condition (3.3a) can be used to simplify (3.4b): 
= => d e( p = d z (p = => d r xp = 



at r = 1 



(3.4b') 



On the (simply-connected) disks, we impose the boundary conditions in the form 
(2.19), i.e. 

= u z = -A h <p } (3.5a) 

= e z -V xu = -A h xp-^d r (r 2 co±) > at z = ±\ (3.5b) 
= d z u z = -d z A h cp J (3.5c) 

These are equivalent to those on the individual components but easier to implement 
since each of (3.5) involves only one of the potentials. The remaining condition (2.19d) 
required on the two circles is insured by (3.4a). 



3.2 Compatibility condition 



We now turn to the compatibility condition (2.14d) for the hydrodynamic problem in 
our potential formulation, where e = e z , n = e r , and is the r = 1 boundary: 

= e r g = e r Vxf = e r Vx ^aj-fe-^ju + s) at r = l (3.6) 

Because e r ■ V x involves only dg and 9 Z , derivatives parallel to the r = 1 boundary, it 
vanishes for all terms in f which are zero or constant at this boundary. For homoge- 
neous boundary conditions (3.4) on the outer cylinder, this is true for 9ju and for S 
defined in (2.6a) in the absence of a magnetic field, leaving only the Laplacian term. 
Referring to (3. Id), we have 

1 

e r • V x Au = d rz Ait d e AA(p (3.7) 

r 

Conditions (3.3a), (3.4c) and (3.5) allow the replacement of Aip and Acp at r = 1 by A^xp 
and Afrcp, which already appear in the governing equations (3.2), leading to 

1 

= d rz A h ip - -d AA h (p at r = 1 (3.8) 
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The complete set of conditions to be imposed on the velocity is then: 



-dgtp + d rz (p = 
r 

d r ip = 
A h( p = 
(p = 

1 

non-axi: 3 rz A^i/7 — -dgAA^cp = 
axi: i/7 = 



> at r = 1 



at r = 



(3.9a) 

(3.9b) 
(3.9c) 
(3.9d) 

(3.9e) 
(3.9f) 



at 



^2 



d z A h cp = 
A h (p = 

These conditions are imposed on (p and xp via the influence matrix method, as will be 
explained in section 4. 



(3.10a) 

(3.10b) 
(3.10c) 



In equations (3.9), we have marked conditions (3.9d) and (3.9e) as applying only to 
axisymmetric or to non-axisymmetric modes. This will be explained in the following 
section. 



3.3 Spatial discretization and symmetry 



We use the spectral spatial discretization: 

L¥J 

/M, Z )= e r^y 

«=-L¥J 

and similarly for (p. The basis functions in the axial direction z are the standard Cheby- 
shev polynomials T^lz/h). Those in the radial direction r are the non-standard poly- 
nomial basis Q„(r) developed by Matsushima and Marcus [22]. Their principal prop- 
erty is that Q„(r) ~ r m as r — >• 0, insuring their regularity at the origin. The basis func- 
tions in the azimuthal direction 6 are the Fourier modes e imB . In (3.11), we do not intro- 
duce new notation for Fourier coefficients, or for coefficients in the 3D tensor-product 
basis, instead distinguishing between physical space values and spectral space coeffi- 
cients by the number and type of superscripts and subscripts. 

The decomposition (3.11) leads to problems and boundary conditions which are de- 
coupled for each Fourier wavenumber m. In fact, because of the reflection symmetry in 



LtJ K-l 2N-1 fry N 

= E E E fknQn(r)T k f )e imd (3.11) 

m=- I ¥ I k=0 n=\m\ \ n / 

L J n+m even 
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z, the problems can be further reduced. A vector field is reflection-symmetric in z if u r , 
Uq are even in z and u z is odd in z, i.e. if the potential ip is even in z and the potential 
<p is odd in z, as can be seen from (3.1). We denote these functions as having parity 
p = s. Quantities related to anti-reflection-symmetric vector fields, i.e. with ip odd and 
<p even, will be denoted as having parity p = a. The boundary conditions (3.9) at r = 1 
can be considered as applying separately to fields of each parity; note that (3.9a) and 
(3.9e) couple potentials of the same parity. The conditions (3.10) at z = ±h/2 can be 
reformulated to apply only to fields of a single parity; for example Ah<P( z = ±|) = 
can be rewritten as A/,^(z = h/2) ± Af l cp(z = —h/2) = 0. Essentiallly, a problem posed 
over the entire cylinder can be viewed as 2(M + 1) problems in the two-dimensional 
meridional half -slice <r < 1,0 <z <h/2. 



Special conditions are applied to the axisymmetric modes. The gauge freedom (2.10b) 
for ip requires the specification of a single value of xp at each z. In (3.9f), we have chosen 
to specify this value at the origin: 

ip(r = 0,9,z) = J^ip m (r = 0,z) e imd = tp m=0 (0,z) (3.12) 

m 

Condition (3.9f) is applied only to the axisymmetric mode, since only this mode con- 
tributes to the sum (3.12). 



For the axisymmetric modes, two important consequences are derived from the calcu- 
lation for an arbitrary function / m=0 (r) 

(rd r f°)(r = R) = rd r f = [* dr d r rd r f° = [* r dr -d r rd r f° = f \ dr A h f° (3.13) 

r=0 Jo Jo r Jo 

Going from left to right in (3.13), one obtains the classic solvability condition required 
by Neumann boundary conditions, since setting the value of d r f°(r = R) is equivalent 
to an integral constraint on A/,/ . In particular, the Neumann boundary condition (3.9b) 
on ip must be replaced by the integral constraint like (3.13) for the axisymmetric mode. 
Going from right to left in (3.13) leads to the conclusion that the only axisymmetric har- 
monic function on a disk that includes the origin is a constant, since A^/° = over [0, R] 
implies d r f°(r = R) = for each R. This implies that the Neumann boundary condition 
in (2.17) is unnecessary to guarantee g = for the axisymmetric mode, and hence that 
the compatibility condition (3.9e) should not be imposed on the axisymmetric mode. 



Corresponding to the removal of the compatibility condition, Marques [1] showed that 
the system of equations governing the axisymmetric modes is of lower order. The cal- 
culation 

- \rf = => f = f = (3.14) 

demonstrates the invertibility of 9+, or equivalently, the impossibility on a disk of a 
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non-zero divergence-free axisymmetric radial vector field which is regular at the ori- 
gin. Using (3.14), equations (2.4)-(2.5) become PDEs of lower order in d r ip° and d r (p°: 

(d t -Re- l A+)d r f u = e e -S° u (3.15a) 

(d t - Re" 1 A+) A+d r <p° u = -e e • (V x S°) (3.15b) 

(d t - Rm- 1 A+) d r (p° B = e e -S° B (3.15c) 

(d t ~ Rm _1 A + ) d^g = e r VxS| (3.15d) 

where A + = 9 r 9+ + d\. In the interests of uniformity we continue to solve the same 
equations for the axisymmetric as for the non-axisymmetric modes, altering only the 
boundary conditions. 

4 Nested Helmholtz and Poisson solvers 

4.1 Temporal discretization 

We briefly mention some aspects of our temporal discretization. A more extensive de- 
scription of both the temporal and the spatial discretization is given in [21, 23]. We 
recall the equations governing the velocity potentials: 

(3* ~ Re~ 1 A)A h tp = e z • V x S =S^, (4.1a) 

(d t - Re~ 1 A)AA h (p = -e z • V x V x S = (4.1b) 

Evolution equations such as (4.1) are typically discretized in time via an implicit scheme 
for the diffusive terms and an explicit scheme for the nonlinear terms. For example, 
with the simplest choice of the backwards and forwards first-order Euler formulas, the 
diffusion equation 

(d t - Ke _1 A) f = S (4.2a) 

becomes 

m±ipm_ Re - lAf{t+st)=s{t) 

(i - £a) /(( + St) =/(()+ St S(f) (4.2b) 



Thus implicit-diffusive /explicit-nonlinear temporal discretization transforms the pa- 
rabolic equation (4.2a) into the Helmholtz problem (4.2b) for f(t + St). Similarly, the 
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temporally discretized versions of the more complicated equations (4.1) give ip(t + St) 
and <p(t + St) as solutions to a sequence of nested Helmholtz and Poisson problems. 
We will not distinguish between the continuous-time parabolic operators of type (4.2a) 
and the discretized Helmholtz operators of type (4.2b) and refer to both as Helmholtz 
problems. 

4.2 Substitution ofDirichlet boundary conditions 

Equations (3.9)-(3.10) give the set of boundary conditions which is to be imposed on 
(4.1). The major difficulty of the poloidal-toroidal formulation is that set (3.9)-(3.10), 
while appropriate for the entire problem, does not provide separate boundary con- 
ditions appropriate to each individual Helmholtz and Poisson problem. Some of the 
conditions involve both if> and (p. Even conditions involving only one potential can 
be problematic because the order of the equations and of the boundary conditions do 
not match. The prototypical example of this occurs in the 2D streamfunction-vorticity 
formulation. At each timestep, one would like to solve successively the Helmholtz 
problem for the vorticity and the Poisson problem for the streamfunction. However, 
no boundary conditions are available for the vorticity, while the streamfunction must 
satisfy both Dirichlet and Neumann boundary conditions. 

The influence matrix technique [3] calls for replacing the problematic boundary con- 
ditions by conditions which are easier to implement numerically, in this case Dirichlet 
boundary conditions on a set of intermediate fields. The values used in these bound- 
ary conditions are determined in such a way that the exact boundary conditions are 
satisfied. We show below the sequence of problems with their associated boundary 
conditions: 




(4.3a) 



U = -- d r(r 2 co±) 



at z = ±- 
2 

at r = 1 



(4.3b) 



axi : 




(4.3c) 



nonaxi : 



fa = ^/(z) 



at r = 1 




(4.3d) 



nonaxi : 



axi : 



AfcV = fa 
i/7 = 

a r i/7 = o 



at r = 
at r = 1 



(4.3e) 

(4.3f) 
(4-3g) 
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(d t - Ke _1 A) g(p = S<p 

4{r)=dzU\ z=± h=0 



at z = ±- 
2 



g<P = W) 



at r = 1 



-dgip + d rz (p 



= 



(4.3h) 
(4.3i) 

(4.3j) 



r=l 



= 



at r = 1 

, h 

at z = ±- 
2 



at r = 1 



(4.3k) 
(4.3£) 
(4.3m) 

(4.3n) 
(4.3o) 



We have introduced intermediate variables ^ and ftp, and required them to obey 
Dirichlet boundary conditions with unknown values oy(z), o~ g (z) and c^{r), or 0. We 
have also introduced the notation Cf(z), c g (z), and c^(r) for quantities which should 
be zero if the actual boundary conditions were satisfied. The boundary conditions in 
(4.3) are identical to (3.9)-(3.10), restated where possible in terms of fy, and The 
influence matrix establishes the correspondence between {oy,c^,(r^} and {cf,Cg,c^}. 
No significance should be attached to the choice of equation in (4.3) at which each c is 
defined, i.e. the elliptic problem with which each of the original boundary conditions 
has been associated. The correspondence serves merely to establish that the number of 
unknown Dirichlet values a is the same as the number of boundary conditions c = 0. 



In order to simplify the notation, we have suppressed the indices labelling the az- 
imuthal Fourier wavenumber m and axial parity p G {s,a}. Each equation in (4.3) 
should in fact be interpreted as applying separately to modes with different (m,p) 
values. Wherever it occurs, dg should be interpreted as multiplication by im, while in 
equation (4.3b), the right-hand-side is axisymmetric and hence should be interpeted 
as zero for m/0. Note that the boundary conditions for the axisymmetric and non- 
axisymmetric modes differ slightly, as explained in section 3.3. 
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4.3 Influence matrix method 



System (4.3) is solved by generalizing the standard decomposition of a linear bound- 
ary value problem into particular and homogeneous problems, in which the bound- 
ary conditions or right-hand-side are set to zero, respectively. Here, the nature of the 
boundary conditions and the intermediate solutions to which they are applied are also 
changed. Historically, the name capacitance matrix has also been used to denote what 
we call the influence matrix; this has guided our choice of notation C. The steps for 
carrying out the influence matrix technique are as follows. 

Preprocessing step (homogeneous solutions): 

• We calculate solutions to the homogeneous problem (S = 0) with a complete set 
of Dirichlet boundary conditions corresponding to the spectral discretization (3.11). 
Specifically, for each Fourier mode m G {0,...,^} and axial parity p € {s,a}, the 
boundary values {cy(z), Cg(z), or^(r)} are set successively to: 

{<r f (z)=T k (%), 
{<7(z)=0, 

{<7(z)=o, 

{cr / (z)=0, 

• For each homogeneous solution, the values of the unsatisfied conditions c p c g , 
are calculated on the boundary. 

• These are collected to form the 2(M + 1) influence matrices C mv , each of size (K + 
N) x (K + N). Each set of Dirichlet boundary values leads to one column of C mp . 

• The influence matrices are inverted to form (C mp ) _1 . Difficulties and techniques re- 
lated to this inversion are discussed in appendix A. 

Each timestep (particular and final solutions): 

• We calculate the particular solution, i.e. the solution to the inhomogeneous problem 
(S 7^ 0) with homogeneous Dirichlet boundary conditions ay = o~ g — — 0. 

• We calculate the values of the unsatisfied conditions {cf,Cg,c^} on the boundary. 
These are separated according to Fourier mode m and axial parity p. 

• Each set (m,p) of c values is multiplied by the corresponding matrix (C m P)~ l to 
obtain appropriate values of {oy,c^,(r^}. 

• The inhomogeneous problem is then solved again with corrected inhomogeneous 
Dirichlet boundary values. 



<r g {z)=0, c7+(r)=0, ^(r)=0| (4.4a) 

o- § (z) = T k (%) , a+(r)=0, <r-(r)=o} (4.4b) 

r g {z)=0, o-+(r) = Q™(r), % {r) = <%{r)} (4.4c) 

0, <r+(r) = QX(r), cr~(r) = -Q™(r)} (4.4d) 



(To Z 
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Axial symmetry is taken into account by examining (4.3). For solutions which are 
reflection-symmetric in z (p = s), Hf is even in z, i.e. k takes only even values in (4.4a), 
while Cg is odd in z, so k takes only odd values in (4.4b). Additionally only (4.4d) is 
used. The corresponding Cf is odd, and c g , c v g m are even in z. The opposite holds for 
fields which are anti-reflection symmetric in z (p = a): k takes only odd values in (4.4a) 
and even values in (4.4b) and only (4.4c) is used. 

This decomposition can be expressed mathematically as a version of the Sherman- 
Morrison-Woodbury formula [3]. Here, a large problem coupling fxp,^,gip,f<p,(p is de- 
coupled by a transformation (the change in boundary conditions) of low rank (K + N). 
The solution to the coupled problem can be obtained from that of the decoupled prob- 
lem using an additional multiplication by a matrix of dimension K + N. 



5 Towards an MHD solver 

We now address the solution of the induction equation in a finite cylinder of finite con- 
ductivity surrounded by a vacuum extending to infinity. For a sphere or axially infinite 
cylinder, the boundary between interior and exterior domains is associated only with 
the radial coordinate. The boundary surrounding a finite cylinder, however, is speci- 
fied as a relation between r and z. One approach is to define the induction equation in 
an integral formulation. The most important advantage is that no boundary conditions 
must be specified. Using this formulation [24], a stationary kinematic dynamo problem 
was solved in a cylindrical geometry. In [25, 26], a finite volume method is used to dis- 
cretize the solution in the interior, which is matched to that in the exterior vacuum via a 
boundary element method. An integral equation formulation was applied to the entire 
domain in [27], and [28] uses finite elements with a penalty method to apply boundary 
conditions. To our knowledge, however, there exists as yet no method applicable to the 
spectral formulation in a finite cylinder. 

5.1 Matching conditions and gauge 

In the remainder of this section, fields or potentials without subscripts or superscripts 
will be taken to refer to the interior magnetic field, while fields or potentials relating to 
the field in the exterior vacuum will be designated by a superscript, e.g. B vac , <p vac . We 
recall the equations governing the interior magnetic potentials: 

(3 f - Rm- 1 A)A h (p = e S = S,;, (5.1a) 
(d t - Rm- 1 A)A h tp = e • V x S = (5.1b) 
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System (5.1) requires two boundary conditions at each bounding surface and supple- 
mentary gauge and compatibility conditions at the horizontal boundary The exterior 
magnetic field in a vacuum is described by a single harmonic potential which requires 
one boundary condition at each bounding surface; see 5.2. Thus a total of three match- 
ing conditions must be applied at each bounding surface: 



o = (B r - K c ) = -M + MM - f ac ) 

O=(B - B v e ac ) = -d r1 p + -d e (d z <p - <p™ 
= (B z - B v z ac ) = -A<p + d z (d z <p - f ac ) 



> on 9Q 



(5.2a) 

(5.2b) 
(5.2c) 



On the bounding cylinder r = 1, equations (5.2b) and (5.2c) (but not (5.2a)) can be 
simplified by choosing the gauge: 



= (d z< p - (p mc ) at r = l 



(5.3) 



leading to: 



= -d e xp + d r (d z (p - cp l 

= d r ip 
= A(p 



at r = 1 



(5.4a) 

(5.4b) 
(5.4c) 



On the disks z = ±h/2, we use the boundary condition (2.19b), which becomes: 

= e 2 • V x (B — B vac ) = -A h ip at z = ±| (5.5a) 

since the exterior magnetic field is curl-free. Thus tp is harmonic on the disks, with 
homogeneous Neumann boundary condition (5.4b) and is therefore constant on each 
disk. The matching conditions (5.2a) and (5.2b) can be applied at the disks to show 
that d z (p — <p vac is constant on each disk, while the gauge condition (5.3) shows that the 
constant is zero: 

= d z (p-(p mc at z = ±^ (5.5b) 
The matching conditions at z = ±| are completed by applying (5.2c) at the disks: 

0= -A(p + d z (d z (p-(p mc ) = -A h (p-d z (p mc at z = ±| (5.5c) 
The final set of gauge and matching conditions to be imposed is (5.3), (5.4) and (5.5). 
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5.2 Exterior magnetic field 



Since the exterior magnetic field has both zero curl and zero divergence, and the ex- 
terior domain is simply connected, (2.11) states that the exterior magnetic field can be 
represented as: 

B mc = V(p vac (5.6) 

where <p vac satisfies: 

A(p mc = outside the cylinder (5.7a) 

Vf ac = at infinity (5.7b) 

Equation (5.7b) supplies the boundary condition on (p vac at infinity, while conditions 
on the cylindrical boundary are provided by coupling with the interior potentials via 
the gauge and matching conditions (5.3), (5.4a), (5.5b) and (5.5c). 

Although (5.7) is posed in the infinite domain outside the cylinder, we can avoid dis- 
cretizing the infinite domain and solving numerically by using known analytic solu- 
tions to the Laplace equation. Our approach is to formulate a complete set of analytic 
solutions to (5.7), each of whose derivatives can be calculated. The exterior solution 
(p vac can be expanded in this set, with coefficients related to values on the cylindri- 
cal boundary. The normal derivatives at the boundary can then be evaluated in terms 
of these coefficients. This defines a correspondence between a set of boundary val- 
ues {<p vac \da\ and a set of normal derivatives {n • V(p mc \sci}. This correspondence, 
or influence matrix, constitutes a basis for the Dirichlet-to-Neumann mapping for the 
domain outside the finite cylinder. The normal derivatives appearing in the match- 
ing conditions (5.4a) and (5.5c) can then be replaced by functions jF r = d r cp vac \ r= i and 
= dz(p vac \ z= ±h of the boundary values {(p vac \dci}- Equations (5.3) and (5.5b) in turn 
relate the boundary values of the exterior and interior potentials via <p vac \dn — ^z^lan- 
The exterior magnetic field no longer appears and the interior problem is closed. Es- 
sentially, we seek to replace the matching conditions (5.3),(5.5b),(5.4a), (5.5c) by: 

0=-d e tp + d rz cp-f r ({d z (p\ m }) at r = l (5.8a) 

r 

= -A h( p-F±({d z( p\ m }) at z = ±\ (5.8b) 

The task is now to obtain a well conditioned matrix representation of the Dirichlet-to- 
Neumann mappings jF r , J 7 ^. 

We have considered two sets of solutions to (5.7). The first set is constructed from the 
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classic spherical harmonics. That is, we expand cp mc as: 

f ac = f»c + £ £ & -C+l)p Zm (cosO^ 

W (5 9) 

p = vV + z 2 , £ = tan" 1 Q 

where P/ m are the associated Legendre polynomials. According to (5.6), <p vac is defined 
only up to a constant, which we may choose such as to set <p™ c = 0. (In two dimensions, 
i.e. for a function whose gradient decays as the cylindrical radius r, rather than the 
spherical radius p, tends to infinity, logarithmic functions would have to be included 
in the expansion because it cannot be assumed that (p vac tends to a constant at infinity.) 
For a given longitudinal Fourier mode m, the solution (5.9) has degrees of freedom 
associated with index I, associated with the latitude. In the standard spectral-physical 
space duality, the set of coefficients (pf™ corresponds to the set of values §™ c {ri,Z[) for 
{ji,Z[) on the boundary and can be determined from them by solving: 

|m|+L-l 

<PT(n,Zi)= E CV (m) ^Wcos6) 

l= H (5.10) 



where L is the number of points on the boundary. Expansion (5.9) readily yields the 
normal derivatives 9 r ^^f c (r ; -,Z;),9 2 ^^f c (r;,Zf) at the boundary in terms of the coeffi- 
cients (pi£f. However, this approach is not feasible, in part because the transform (5.10) 
is extremely poorly conditioned, like all other transforms involving monomials. As I 
increases, the functions p\ become spiked at the largest values of p\ and zero else- 
where, a difficulty which does not arise on a spherical surface where p is constant. 

We have also considered a second set of solutions to (5.7), constructed from the equally 
classic free-space or fundamental Green's functions: 

where cr(x') is a distribution on the cylindrical surface 9Q which is calculated in such a 
way as to yield a particular set of boundary values for <p vac (x) and G(x,x') = 4zr/ |x — 
x'| satisfies 

A' x G(x,x') = S(x-x r ) (5.12) 

The derivatives of expansion (5.11) are obtained by differentiating the Green's func- 
tions: 

V^xH/^x'j^^' M (5.13) 
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This approach is discussed and shown to perform quite well for a two-dimensional 
test problem in [29]. The exterior fields generated approximate the solution uniformly 
near the boundary and converge exponentially with resolution; the influence matrix 
representing the Dirichlet-to-Neumann mapping is well conditioned. 



5.3 Magnetic compatibility condition 

The magnetic compatibility condition, obtained by substituting (2.12c) into (2.14d), is: 

= e r ■ g B = e r • ((3, - Rm- l ^j B + S B ) at r = l (5.14) 

Contrary to the velocity, the magnetic field does not vanish on the boundary so 9f B 7^ 0. 
The nonlinear term Sg = V x (u x B) does vanish at the boundary since u| r= i = and 
its radial curl e> • Sg contains no normal derivatives of u. The remaining terms are 
evaluated using (3.1a) and (3.1c): 

= e r -(d t -Rm- 1 A)B = ^d e (d t -Rm- 1 A)ip + d rz (dt-Rm- 1 A)(p at r = l (5.15) 

The magnetic compatibility equation must use the same time discretization as the evo- 
lution equations, here backwards Euler. Although the the time derivative in (5.15) may 
seem difficult to include in an implementation of the influence matrix, the boundary 
operator in (5.15) can be decomposed into two parts, one which acts on the homoge- 
neous solution at time t + St (and contributes to the influence matrix) and the other 
which acts on the particular solution at time t + St and the actual solution at time t. 
Details are given in [21]. 

The velocity compatibility condition (3.9e) must also be modifed to include a contribu- 
tion from the magnetic field: 

1 

= d rz A h ip u - -d 9 AA h (p u — e r • V x (B • V)B at r = 1 (5.16) 



5.4 Nested elliptic problems and influence matrix 

The set of nested Helmholtz and Poisson problems (5.1), together with the conditions 
(5.4b)-(5.4c), (5.5a), (5.8) and (5.15) can be solved using the influence matrix technique, 
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as was done for the velocity in section 4.2 and in equation (4.3). 



(d t - Re- 1 A") f = 



f = 

axi : r drf = 
Jo 

nonaxi : / = crAz) 

it 

c f (z) = 







, h 

at z = ±2 
at r = 1 

at r = 1 



-d e (dt - Rra _1 A)i/> + d rz {d t - Rm _1 A)0 



(5.17a) 

(5.17b) 
(5.17c) 

(5.17d) 







r=l 



A h xp = f 

axi : i/7 = 
nonaxi : d r ip = 



at r 
at r 




1 



(5.17e) 
(5.17f) 
(5.17g) 



(ai-Re- 1 A)g = S (f 

ft 

Cg(z) =Atf>| r= i=0 
ft 



at r = 1 



at z = ± - 



4 (r) = -g| z=±S " J?({3z^|an}) = 



(5.17h) 
(5.17i) 

(5.17J) 



&h<P = g 
<p = o-<p(z) 

it 

c<p(z) ; 



-a i/7 + a rz </> 



at r = 1 



(5.17k) 
(5.17*) 



r=l 



The major differences between the velocity and magnetic cases are the reduction from 
five to four in the number of elliptic problems, the time derivative in the magnetic com- 
patibility equation, and the presence of the Dirichlet-to-Neumann mappings !F r and 
Tf- which have replaced the normal derivatives of the exterior solution <p vac . We recall 
the meaning of F r {{d z <p\dri\) an d ^t({^z(p\dci})'- the set {^z^lan} provides Dirich- 
let boundary values for the exterior Laplace problem and T r and are the normal 
derivatives of the exterior solution at the boundaries. 

The Dirichlet-to-Neumann mappings, as well as each solution and problem in (5.17), 
decouple according to azimuthal Fourier mode m and axial parity p £ {s,a}. As before, 
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we construct the influence matrix by solving homogeneous versions of (5.17), with 
and Sq set to zero in (5.17a) and (5.17h) and a complete set of Dirichlet boundary values 
Up <jg, a^, (T<p. Evaluating Cf, c g , c^, c<p yields the influence matrix. 



6 Conclusion 



We have presented a poloidal-toroidal formulation for solving the time-dependent 
three-dimensional magnetohydrodynamic equations in a finite cylinder. While pre- 
serving the original mathematical formulation described in [1] and later tested on a 
linear stability Rayleigh-Benard convection problem in [2], we incorporated the influ- 
ence matrix technique [3] for decoupling the boundary and compatibility conditions 
emerging from the potential formulation. We have also described an extension of this 
algorithm to the induction equation governing the evolution of the magnetic field. 

The most important advantage of using the toroidal-poloidal decomposition is that 
the divergence-free character of the velocity and magnetic fields is imposed exactly, by 
construction. For the induction equation, the potential formulation makes it possible to 
solve for the magnetic field without introducing an artificial numerical magnetic ana- 
logue to the hydrodynamic pressure which has no physical meaning. In addition, us- 
ing scalar functions instead of components of vector fields simplifies and homogenizes 
the usage of differential operators. The influence matrix technique allows the poloidal- 
toroidal formulation to be sufficiently economical to be used for time-integration. 

We have implemented and validated this method for the hydrodynamic von Karman 
problem of flow in a cylinder driven by counter-rotating disks, using a spectral dis- 
cretization which is regular on the cylindrical axis. These results are presented in a 
companion article [23]. In extending this algorithm to the full magnetohydrodynamic 
problem, no difficulty is posed by the induction equation, whose structure is simpler 
than that of the Navier-Stokes equation. Instead, the main difficulty is that the mag- 
netic field is not specified at the domain boundary but must instead satisfy matching 
conditions between the interior domain (here a finite cylinder) and the exterior domain 
(here an infinite vacuum). We have developed a formalism involving the Dirichlet- to- 
Neumann mapping for eliminating the exterior magnetic field, which has been imple- 
mented and validated for a two-dimensional test problem [29]. Future research will 
focus on implementing the poloidal-toroidal formulation for the full magnetohydro- 
dynamic problem. 
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A Regularizing the influence matrix 



In section 4.3, we have assumed that the influence matrices C vm are invertible, which is 
actually not the case. The influence matrices are non-invertible for several reasons. The 
first issue is geometric. The finite cylinder has corners at which conditions at z = ±| 
and r = 1 must both be satisfied. When formulated in spectral rather than physical 
space, the redundant conditions correspond to a linear combination of rows and cannot 
be easily identified. A second factor is the discretization of the Poisson and Helmholtz 
solvers, in particular the replacement of the highest-wavenumber equations by the 
boundary conditions mandated by the r method. A third cause is the decrease in poly- 
nomial order due to differentiation by boundary operators. For numerical reasons, the 
eigenvalues corresponding to these directions may be nearly zero, rather than exactly 
so. 

One remedy [3] consists of thresholding: diagonalizing C vm and replacing the eigen- 
values whose absolute values are below a certain experimentally-determined thresh- 
old by an arbitrary value, say 1, leading to an invertible matrix. The justification of 
this manipulation of the spectrum is that the eigenvectors corresponding to the zero 
eigenvalues play no role in satisfaction of the boundary conditions. This is true if the 
linear system of equations defined by the influence matrix and the right-hand-side is 
underdetermined, i.e. if the right-hand-side belongs to the image space of the influence 
matrix. Because the particular solutions are determined using the same nested solver 
used for constructing the homogeneous solutions, this is in fact the case. 

However, a major problem remains. Even after eigenvalues are eliminated which would 
be exactly zero if infinite precision were used, the resulting matrices still have very 
small eigenvalues, i.e. they are still poorly conditioned. There are various causes for 
this. Some boundary value distributions are almost linearly dependent. More impor- 
tantly, because some boundary conditions are of higher differential order then others, 
the magnitudes of different portions of the influence matrices are very different. We 
shall call these eigenvalues small, in contrast to those which would be zero in infinite 
precision, which we shall call simply zero eigenvalues. The small eigenvalues depend 
on the spatial resolution and on the product Re /St, a parameter which appears in the 
Helmholtz problem. As the resolution or Re /St are increased, an increasing number of 
small eigenvalues appear, whereas the number of zero eigenvalues depends only on 
the geometry and on the kind of boundary conditions. 

The condition number, approximately the ratio between the largest and smallest eigen- 
values values of a matrix, is an upper bound on the number of lost meaningful digits in 
the numerical solution to a linear equation involving this matrix. If the threshold is 
chosen such as to lower the condition number to an acceptable value (O(10 8 ) — O(10 10 ), 
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then small eigenvalues are eliminated, in addition to zero eigenvalues, leading to errors 
in satisfaction of the boundary conditions that are much higher than machine precision. 
The challenge is, first, to distinguish between the zero and small eigenvalues so as to 
eliminate only the zero eigenvalues and, second, to improve the condition number of 
the adjusted matrix in some way other than by eliminating the small eigenvalues. 

We first modified the thresholding procedure by using the singular value decompo- 
sition (SVD) rather than diagonalization. The advantages of the SVD is, first, that it 
always exists and, second, that the matrix of singular vectors is better conditioned 
than the eigenvector matrix, because the left and right singular vectors are orthogonal, 
in contrast to eigenvectors, which may be close to linearly dependent. As an example, 
we consider the influence matrix for spatial resolution (N = 96) x (K = 192), Reynolds 
number Re = 10 4 and time step 5t = 10~ 2 . The magnitudes of the singular values 7/ of 
the influence matrix block C 1,s with azimuthal Fourier wavenumber m = 1 and axial 
parity p = s are presented on figure 3a. C 1,s has one zero singular value (i.e. a singular 
value which would be exactly zero in infinite precision). Figure 3(left) shows that this 
value is numerically 10~ 21 , separated from the next smallest singular value. In contrast, 
the zero eigenvalue and smallest remaining eigenvalue are of the same size and hence 
cannot be distinguished. However, thresholding, whether by replacing this singular 
value or the corresponding eigenvalue, still leaves the condition number unacceptably 
high, on the order of 10 6 /10~ 15 = 10 20 for this example. 




50 100 150 200 250 50 100 150 200 250 

i i 

Fig. 3. Singular values {7,} for case with (N = 96) x (K = 192), Re/5t = 10 6 , m = 1, p = s. Left: 
original influence matrix C 1,s . Right: scaled influence matrix (C 1,s )' = [ a ] C 1,s [ jS ]. 



The matrix condition number can be decreased more effectively by scaling prior to 
thresholding. If each row is divided by its norm, the condition number of the matrix 
is significantly reduced, down to 10 n f or m > and below 10 8 for m = 0. Additional 
scaling of columns does not significantly change the condition number. The condition 
number of the m = matrices are sufficiently decreased by scaling only their rows, be- 
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cause their poloidal and toroidal potentials are not coupled. For m > 0, the influence 
matrices can be further improved by scaling blocks corresponding to different combi- 
nations of types of test functions and boundary conditions. The influence matrices C pm 
f or m > are composed of 9 submatrices: the three columns correspond to the three 
types of imposed simplified Dirichlet boundary conditions {ag(z) l af(z) l a g h (r)} while 
the three rows correspond to the three quantities {c g (z),Cf(z),c g L (r)} reflecting the ac- 
tual boundary conditions. We represent the | • |oo norm of the corresponding submatrix 
of C? m by Cjj, where, for example, C\j_ designates the norm of the (c^,ay) submatrix. 

For the resolution (N = 96) x (K — 192), the matrix norms cu of each of the blocks are: 
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C23 
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C31 
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We now wish to scale the block-rows and block-columns in such a way as to make the 
norms c'- { of the resulting scaled blocks equal to one another. 
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In general there exist no 
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= c'22 



(A.2) 



c 32 = c 33 = 1. We can instead require 



'23 



0L\$\C\\ = CtlfilCll = «3jS3C33 = 1 
Ct2plC21 = 0i\^>2C\2 
«3j6lC3i = 0C^ 3 C 13 



(A3) 



The system (A.3) has an infinite number of possible solutions, from which we can select 
the following: 



Oil = 

h = 



C21C32C33 
C11C12C23 

C11C23 
C11C21C32C33 
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C32C33 
C22C23 
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C 23 R _ 1 
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C22C32C33 C33 



(A.4) 
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After scaling using (A.2)- (A.4), the influence matrix C 1,s has the following structure: 



c n ; c i2 ; c i3 



1 ;10 



2 



10 



-2 




c 21 ! c 22 ! c 23 



10 



1 



1 



(A.5) 



c 31 ; c 32 ; c 33 







1 



1 



The zero singular value is then easily identified and replaced by 1, leading to a condi- 
tion number of 10 11 , like that for simple row or column scaling. But if the block scaling 
is followed by row scaling and then replacement of the zero singular value, then the 
condition number is further reduced to the acceptable value of 10 8 . The singular values 
of the matrix after block and row scaling are presented on figure 3(right). 

We note that operators with high condition numbers are inherent in the numerical 
discretization of partial differential equations; for example, the ID second derivative 
operator with homogeneous boundary conditions using a basis of K Chebyshev poly- 
nomials and corresponding grid has condition number 0(K 4 ). The requirements for 
the solution of the linear systems that occur in this context are not those of numerical 
linear algebra: the right-hand-side is not arbitrary, but results from time-integration, 
and not all components are of equal weight. In our case, we find that after an initial 
small integration time of T = 100 St, the right-hand-side is always such that the influ- 
ence matrix, scaled and regularized to reduce its condition number to O(10 8 ), can be 
inverted to satisfy the constraints in system (4.3) to machine accuracy [23]. 

One welcome consequence of scaling is that it separates the zero singular values which 
result from non-invertibility of the matrix from the small singular values which result 
from poor conditioning of the various components of the influence matrix). Without 
scaling, the number of singular values below a fixed threshold depends on the spatial 
resolution and so the zero eigenvalues cannot be reliably identified and removed. More 
importantly, scaling vastly improves the condition number of the influence matrix, 
insuring satisfaction of the constraints to machine accuracy. 

References 

[1] F. MARQUES, On boundary conditions for velocity potentials in confined flows: 
Application to Couette flow, Phys. Fluids A 2, 729 (1990). 

[2] F. Marques, M. Net, J. M. Massaguer & I. Mercader, Thermal con- 
vection in vertical cylinders. A method based on potentials, Comput. Meth- 
ods Appl. Mech. Engrg. 110, 157-169 (1993). 



[3] L. S. TUCKERMAN, Divergence-free velocity fields in nonperiodic geometries, 
J. Comput. Phys. 80, 403-441 (1989). 



28 



[4] D. REMPFER, On boundary conditions for incompressible Navier-Stokes prob- 
lems, Applied Mechanics Reviews, 59, 107-125 (2006). 

[5] J.M. LOPEZ, F. MARQUES, J. SHEN, An efficient spectral-projection method for 
the Navier-Stokes equations in cylindrical geometries - Three-dimensional cases, 
J. Comput. Phys. 176, 384-401 (2002). 

[6] J.U. BRACKBILL & D.C. BARNES, The effect of nonzero V • B on the numerical 
solution of the magnetohydrodynamic equations, J. Comput. Phys. 35 426 (1980). 

[7] K. Chan, K. Zhang, J. Zou, G. Schubert, A non-linear 3D spherical a 2 dy- 
namo using a finite element method, Phys. Earth Planet. Int. , 128 35-50 (2001). 

[8] P. S. MARCUS, Effects of truncation in modal representations of thermal convec- 
tion, J. Fluid Mech. 103, 241-255 (1981). 

[9] G. A. GLATZMAIER, Numerical simulation of stellar convective dynamos. 1. The 
model and method. J. Comput. Phys. 55, 461-484 (1984). 

[10] M. DUDLEY & R. JAMES, Time-dependent kinematic dynamos with stationary 
flows, Proc. Roy. Soc. London A 425, 407-429 (1989). 

[11] G. A. GLATZMAIER AND P. H. ROBERTS, A three-dimensional self-consistent 
computer simulation of a geomagnetic field reversal, Nature 337, 203 (1995). 

[12] A. TlLGNER, A kinematic dynamo with a small scale velocity field, Phys. Rev. A 
226, 75-79 (1997). 

[13] R. HOLLERBACH, A spectral solution of the magneto-convection equations in 

spherical geometry, Int. J. Num. Meth. Fluids 32, 773 (2000). 
[14] H.B. SQUIRE, On the stability for three-dimensional disturbances of viscous flow 

between parallel walls, Proc. Roy. Soc. Lond. A 142, 621 (1933). 
[15] P. J. SCHMID & D. S. HENNINGSON, Stability and Transition in Shear Flows, 

Springer, 2001. 

[16] J. Antonijoan, F. Marques & J. Sanchez, Non-linear spirals in the Taylor- 

Couette problem, Phys. Fluids 10, 829 (1998). 
[17] A. P. WILLIS & C. F. BARENGHI, Hydromagnetic Taylor-Couette flow: numerical 

formulation and comparison with experiment, J. Fluid Mech. 463, 361-475 (2002). 
[18] T. VON K ARM AN, Uber laminare und turbulente Reibung, Z. Angew. Math. Mech. 

1, 233-252 (1921). 

[19] M. Bourgoin, L. Marie, F. Petrelis, C. Gasquet, A. Guigon, J. B. Luciani, 

M. MULIN, F. NAMER, J. BURGETE, A. CHIFFAUDEL, F. DAVIAUD, S. FAUVE, 
P. ODIER & J. F. PlNTON, Magnetohydrodynamics measurements in the von Kar- 
man sodium experiment, Phys. Fluids 14, 3046 (2002).. 

[20] C. NORE, L. S. TUCKERMAN, O. Daube & S. XlN, The 1:2 mode interaction in ex- 
actly counter-rotating von Karman swirling flow, J. Fluid Mech. 477, 51-88 (2003). 

[21] P. BORONSKI, A method based on poloidal-toroidal potentials applied to the von 
Karman Flow in a finite cylinder geometry, Ph.D. Thesis, Ecole Polytechnique, 
2005. 

[22] T. MATSUSHIMA & P. S. MARCUS, A spectral method for polar coordinates, 



29 



J. Comput. Phys. 120, 365 (1995). 
[23] P. BORONSKI & L. S. TUCKERMAN, Poloidal-toroidal decomposition in a finite 
cylinder: II. Discretization, regularization and validation, submitted to J. Com- 
put. Phys. 

[24] W. DOBLER & K. H. RADLER, An integral equation approach to kinematic dy- 
namo models, Geophys. Astrophys. Fluid Dyn. 89, 45-47 (1998). 

[25] A. B. ISKAKOV, S. DESCOMBES & E. DORMY, An integro-differential formula- 
tion for magnetic induction in bounded domains: boundary element-finite vol- 
ume method, J. Comput. Phys. 7, 540-554 (2004). 

[26] A. B. ISKAKOV & E. DORMY, On magnetic boundary conditions for non-spectral 
dynamo simulations, Geophys. Astrophys. Fluid Dyn. 99, 481-492, (2005). 

[27] M. Xu, F. STEFANI & G. GERBETH, The integral equation method for a steady 
kinematic dynamo problem, J. Comput. Phys. 196, 102-125 (2004). 

[28] J.L. GUERMOND, R. LAGUERRE, J. LEORAT & C. NORE, An interior penalty 
Galerkin method for the MHD equations in heterogeneous domains, J. Com- 
put. Phys. 221, 349-369 (2007). 

[29] P. BORONSKI, Spectral method for matching exterior and interior elliptic prob- 
lems, J. Comput. Phys., in press. 

[30] M.F.M. SPEETJENS, Three-dimensional chaotic advection in a cylindrical domain. 
Ph.D. Thesis, Fluid Dynamics Laboratory - Eindhoven University of Technology, 
2001, ISBN 90-386-1879-4. 



30 



